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I.  INTRODUCTION 


The  purpose  of  this  investigation  has  been  to  improve  a finite 
element  program  which  performs  three  dimensional  analysis  for  impul- 
sively loaded  laminated  plates  (acronym  TIP)1.  This  TIP  program  has  a 
finite  element  model  that  uses  the  quadrilateral  to  define  the  shape 
of  the  element  in  the  plane  of  the  plate  and  then  arranges  a number  of 
these  elements  through  the  thickness  to  describe  the  necessary  number 
of  material  layers.  Each  layer  then  has  its  own  material  properties. 
The  model  is  nonlinear  since  it  allows  for  large  plate  deflections  and 
for  material  yield  effects.  The  original  program  only  allowed  for  iso- 
tropic elastic-perfectly  plastic  solids.  The  dynamic  equations  are 
obtained  by  lumping  the  mass  of  the  plate  into  the  nodal  points  of  the 
finite  element  model  and  then  solving  kinematic  equations  of  motions  by 
use  of  a finite  difference  technique.  The  improvements  to  the  TIP  pro- 
gram were  to  decrease  the  time  that  the  program  ran  by  use  of  large 
time  increments  in  the  finite-difference  equations  and  to  develop  an 
orthotropic  elastic-perfectly  plastic  analysis  which  was  then  further 
improved  to  include  orthotropic  elastic-visco-plastic  analysis. 

The  limitations  on  the  size  of  the  integration  time  step  are  a 
direct  result  of  the  finite-element  model  which  treats  the  plate  as  a 
lumped  mass  system  where  the  individual  masses  are  placed  at  the  nodes. 
Consequently,  as  the  numerical  integration  proceeds  the  masses  move  re- 
lative to  each  other,  and  if  the  integration  time  step  is  too  large, 
an  artificial  oscillation  was  observed  by  examining  the  internal  forces 
acting  on  each  mass  for  successive  time  steps.  When  a relatively  large 
time  step  is  used  these  forces  will  change  sign  for  each  successive 
time  step.  The  worst  condition  occurs  in  the  thickness  direction, 
since  the  masses  will,  in  general,  be  separated  by  a much  smaller  dis- 
tance in  this  direction.  In  addition  to  this,  it  has  been  observed 
from  previous  numerical  results  that  the  two  inplane  displacements  vary 
relatively  smoothly  through  the  thickness  and  therefore,  they  are  modi- 
fied to  vary  linearly  without  restricting  any  sheaT  deformation. 

The  orthotropic  yield  analysis  uses  the  same  technique  that  was 
used  in  the  original  TIP  program  with  only  the  yield  criterion  chang- 
ing. For  the  orthotropic  elastic-perfectly  plastic  analysis.  Hill's 
Yield  Criterion^  is  used  instead  of  von  Mises  Yield  Criterion.  Addi- 
tional orthotropic  yield  stresses  are  also  needed  as  additional  input. 
For  the  elastic-visco-plastic  analysis  a strain  rate  dependence  is 


Zak,  Adam  R.,  "Nonlinear  Dynamic  Analysis  of  Flat  Laminated 
Plates  by  the  Finite  Element  Method,"  Final  Report,  Contract  No. 
DAADO5-73-C-0197,  University  of  Illinois,  February,  1977. 

2Hill,  R.,  "The  Mathematical  Theory  of  Plasticity,"  Clarendon 
Press,  Oxford,  1950. 
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introduced  by  use  of  an  extension  of  the  isotropic  Bingham  material^ 
to  the  orthotropic  case. 

II.  MODIFICATION  FOR  LARGE  TIME  INCREMENTS 


Analysis 

In  this  analysis  the  in-plane  (u  and  v)  displacements  are  handled 
separately  from  the  displacements  through  the  thickness  (w) , but,  both 
utilize  the  finite  difference  equation: 


[<FI}n.l+(?-2)  {V„‘{Vn 

[lFE),..*(r!)fV1*1V,-i] 

where  {A}  is  the  displacement  matrix,  3 is  the  acceleration  parameter, 
h is  the  time  interval,  {Fj}  is  the  internal  force  matrix,  {Fe}  is  the 
external  force  matrix,  [M]  is  the  mass  matrix,  and  the  subscripts  n, 
n-1,  n+1  denotes  time  intervals. 

In  analyzing  the  u and  v displacements,  it  is  assumed  that  they 
are  linearly  dependent  through  the  thickness.  This  forces  plane  sec- 
tions to  remain  plane. 

This  first  assumption  results  in  the  following  equations: 


(A)n+l  - 2{41„  - 
+ 3h2  [M]"1 

+ Bh2  {M]_1 


u = ql  + zq2  (2) 

V " q3  + Zq4 

where  q^,  k=l,  4,  are  unknown  coefficients  called  the  transformed  dis- 
placements and  z is  the  distance  in  the  z-direction  of  the  node  from 
the  center  of  gravity.  The  importance  of  having  z be  the  distance  from 
the  center  of  gravity  will  be  discussed  when  the  transformed  mass  ma- 
trix is  discussed.  In  matrix  notation  equation  (2)  becomes: 


3 

Cristescu,  No.,  "Dynamic  Plasticity,"  North  Holland  Publishing 
Company,  1967. 


6 


{A}  = [TF]  {q} 


(3) 


where  {q}  is  the  matrix  of  transformed  displacements  and  [TF]  is  the 
transformation  matrix  described  below.  Letting  l be  the  number  of 
layers  of  material  and  i=l,  £+1  be  the  nodal  location  in  the  thickness 
direction,  the  transformation  matrix  can  be  written  as: 


[TF]  = 


[TFj] 

[tf2] 

• 

[TF  ] 
• 1 


[TF  ] 
1+1 


where 


[TF.]  = 


1 

0 


z. 

1 

0 


and  Z-  is  the  distance  of  node  i from  the  center  of  gravity.  i=l  is 
the  node  at  the  bottom  of  the  plate  and  i=£+l  is  the  node  at  the  top. 


Since  the  displacements  are  written  in  terms  of  transformed  dis- 
placements, then  the  forces  should  be  written  in  terms  transformed 
forces.  Letting  {fg}  be  the  matrix  of  external  forces  corresponding  to 
{q},  and  {fj}  be  internal  forces  also  corresponding  to  {q},  the  theory 
of  virtual  work  states: 


d{A}T{Fj}  = d{q}T{fj} 


(6) 


Transposing  equation  (3)  yields: 


{A}T  = {q}T  [TF]T 


(7) 


Substituting  equation  (7)  into  (6]  yields: 
d{q}T  [TF]T  {Fj}  = d{q}T  {fj} 
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Therefore: 


{fj}  = [TF]T  {Fj} 

Similarly: 

{fE>  = [TF] T {F£} 


Finding  a transformed  mass  matrix  [m]  starts  with: 


- [M]  {A}  = {F}, 


inertia 


By  virtual  work: 


- d{A}T[M]{A}  = d{A}T{F} 


inertia 


- d{q}T[n,H,}  . d£A>T{F>.nert.a 


d{A)hM]{A}  = d{q}^[m]{q} 


Substituting  equations  (3)  and  (7)  into  (14)  yields: 


d(q)  [TF]T[M][TF]{q}  = d{q}T[m]{q} 


Dividing  out  the  unnecessary  terms  gives: 


[m]  = [TF] 1 [M]  [TF] 


(9) 

(10) 

(11) 

(12) 

(13) 

(14) 

(15) 

(16) 


Now  the  reason  for  Z ^ being  the  distance  from  the  center  of  gravity  will 
become  apparent.  The  calculations  are  much  simplified  by  the  mass  ma- 
trix being  a diagonal  matrix  as  is  the  case  for  the  original  mass  ma- 
trix. The  original  mass  matrix  was: 

[M]  = {M.i  (17) 

where  j varies  over  the  total  degrees  of  freedom. 
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Performing  the  matrix  multiplication  in  equation  (16)  using  equations 
(17),  (4),  and  (5)  produces 


[m]  = 


£+1 

£+1 

EM. 

i 

EM.  Z. 

l l 

0 

0 

i=l 

i=l 

l+l 

l+l 

Jl+1 

EM.  Z. 

l l 

EM.  Z? 

1 X 

EM.Z. 

J T_ 

0 

i=l 

i=l 

i=l 

Jl+1 

l+l 

l+l 

0 

EM.Z. 

l l 

EM. 

l 

EM.Z. 

1 X 

i=l 

i=l 

i=l 

l+l 

l+l 

0 

0 

EM.Z. 
l l 

EM.Z? 

l l 

i=l 

i=l 

But 


£+1 

EM.  Z.  = 0 

l i 

i = 1 


(19) 


by  definition  of  the  center  of  gravity,  thus  causing  [m]  to  be  a diag- 
onal matrix. 


Referring  to  equation  (1)  the  only  element  that  still  must  be 
transformed  is  the  deflections  at  past  time  intervals.  From  equation 
(3): 


{q)n  = [TF]""1{A}n  (20) 

Rather  than  finding  the  inverse  of  the  entire  transformation  matrix,  it 
is  only  necessary  to  get  the  inverse  of  two  of  the  submatrices  in 
equation  (4)  since  it  is  only  necessary  to  solve  for  four  unknowns. 

This  is  easily  done  by  hand  and  results  in: 
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zrz2 


zrz2 


i - 


zrz2 


Vz2 


zrz2 


zrz2 


zrz2 


Zl"Z2 


(21) 


Substituting  the  transformed  quantities  into  equation  (1)  yields: 


- 2 {<>}„  - t^n-l  ‘22> 


f{fT}  * fi-2 

\ {fT}  . + {fT}  01 

L 1 n \B 

} I n-1  I n-2 J 

* 6h2  W*  |-{fE}ntl  * (£  - 2)  «E>n  * 


It  should  be  noted  that  the  internal  transformed  forces  are  displaced 
by  one  time  increment.  Because  these  forces  are  small,  even  at  large 
time  increments,  this  yields  accurate  results.  Equation  (22)  is  thus 
solved  and  equation  (3)  transforms  the  results  to  global  displacements. 

Although  this  time  lag  is  acceptable  for  the  in-plane  displace- 
ments, it  is  not  acceptable  for  the  w displacements.  The  reason  for 
this  is  the  external  force  is  being  applied  in  the  w direction,  thus 
making  these  internal  forces  more  reactive  to  larger  time  intervals . 

In  order  to  account  for  the  change  in  the  internal  force  a model  was 
sought  to  couple  the  deflections  through  the  thickness. 

In  finding  a model  to  represent  what  happens  through  the  thickness, 

it  is  necessary  to  see  what  the  unknowns  are.  From  equation  (1)  the 

unknowns  are  {A}  - and  {FT}  , . All  the  other  terms  are  known.  In 

n+i  I n+1 
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order  to  predict  what  {Fj}  is,  it  is  necessary  to  couple  the  de- 
flections through  the  thickness  and  to  assume  all  strains  small  when 
compared  to  the  strain  in  the  w-direction.  This  can  be  done  by  let- 
ting: 


{FI>„.1  ' fFI}n  * {AFI}n*l 


(23) 


where  (AFj}n+^  is  the  change  of  the  internal  force  between  time  inter- 
vals. The  deflection  in  the  w direction  is  then  coupled  by  the  model 
shown  in  Figure  1.  This  model  assumes  the  stiffness  between  the  nodal 
points  in  the  thickness  direction  is  much  greater  than  the  stiffness 
between  in-plane  nodal  points.  As  long  as  the  external  force  is  in 
the  w direction  this  is  a good  assumption. 


From  Figure  1: 
AF 


(24) 


Ii,n+1  ki  fAi+l,n+l  " Ai,n+P  “ ^Ai+l,n  " Ai,n^j 

" ki_i  [(Ai,n+l  “ Ai-l,n+l5  ' (Ai,n  “ Ai-l,n^] 

where  i refers  to  the  nodal  location  through  the  thickness  and  n refers 
to  the  time  increment.  The  predicted  stiffness  (k^)  is  gotten  from  the 
orthotropic  properties  (C^ , i=l,6,  j=l,6).  In  matrix  notation: 


c -> 

a 

x 


s = 


xy 


xz 


Y z 


C11 

C12 

C13 

0 

0 

0 

< ■> 
e 

X 

C21 

C22 

C23 

0 

0 . 

0 

e I 
y 

C31 

C32 

C33 

0 

0 

0 

< 

e 

z 

> 

0 

0 

0 

C.  . 

44 

0 

. 0 

e 

xy 

0 

0 

0 

0 

C55 

0 

e 

XZ 

0 

0 

0 

0 

0 

e 

66j 

yz 

(25) 


Using  the  assumption  that  all  strains  are  small  when  compared  to  the 
strain  in  the  w-direction  yields 


a 

z 


e 

z 


(26) 
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This  gives  the  stiffness  for  a unit  cube  equal  to  C^* 
Therefore : 


, 33. Area, 

k.  = i i 
1 


(27) 


where  t^  is  the  thickness  of  layer  i.  Area  i is  the  area  used  to  com- 
pute the  mass  of  the  nodal  point,  and  C33i  is  the  orthotropic  property 
of  the  quadrilateral  that  the  node  lies  in. 

Letting: 

a!  . = 2 A.  +A. 

i,h+l  i,n  i,n-l 


(28) 


and  substituting  equations  (23)  and  (28)  into  (1)  produces: 


* 8h2 

Ai,n+1  = Ai,n+1  + M“ 

1 


ki  Ai+l,n+l  Ai,n+1 


(29) 


- k 


i-1 


A.  , - A.  , . - k.  A.  . - A. 

i,n+l  1-1,  n+1  1 i+l, n i,n 


+ ki-l 


A.  - A.  . 

i,n  l-l, n 


where  the  only  unknowns  are  the  deflections  at  time  n+1.  This  pro- 
duces i=£+l  (A  is  the  number  of  layers)  number  of  simultaneous  equa- 
tions which  can  be  solved  for. 


Numerical  Results  and  Conclusions 


The  modified  computer  program  was  applied  to  the  dynamic  plate 
problem  for  which  BRL  experimental  data  is  available,  and  which  was 
analyzed  by  the  original  computer  program. 
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In  the  original  analysis,  before  the  modifications  described  here 
were  incorporated,  the  maximum  time  step  which  could  be  used  without 
numerical  instability  was  0.25  microseconds  (ys) . In  the  modified  pro- 
gram time  step  as  large  as  10  ys  was  used.  Figures  2 to  4 contain  some 
typical  results.  Figure  2 shows  the  results  for  a time  step  of  2.5  ys 
and  Figure  3 has  similar  results  for  5 ys.  Figure  4 contains  the  re- 
sults for  a variable  time  step  calculation  where  1 ys  integration  inter- 
val was  used  for  time  from  0 to  10  ys,  followed  by  a time  step  of  5 ys 
from  10  ys  to  60  ys,  and  finally  10  ys  time  step  for  time  of  60  ys  to 
100  ys.  The  results  show  the  plate  center  deflection  as  a function  of 
time  for  the  bottom  and  the  top  of  the  plate.  These  results  are  com- 
pared with  the  results  from  the  HEMP  solution  and  the  experimental  data 
which  are  only  given  for  the  bottom  of  the  plate.  The  different  time 
steps  were  used  in  order  to  compare  the  effect  of  these  on  the  accuracy. 
It  can  be  observed  that  the  results  for  the  time  step  of  2.5  ys.  Figure 
2,  compare  more  closely  with  the  experimental  data  than  the  results  for 
5 ys.  Figure  3.  The  results  for  a variable  time  step.  Figure  4,  seem 
to  give  even  more  accurate  results  and  compare  quite  well  with  the  HEMP 
calculations  and  the  experimental  data.  It  may  be  mentioned  that  the 
results  for  large  time  steps  introduce  artificial  errors  due  to  large 
distortion  of  the  pressure  load  and,  therefore,  the  errors  may  not  all 
be  due  to  numerical  integration.  However,  the  main  conclusion  is  that 
the  method  has  been  made  quite  stable  and  large  integration  steps  can 
be  used. 


III.  ORTHOTROPIC  ELASTIC-PERFECTLY  PLASTIC  YIELD 


Analysis 

The  orthotropic  analysis  follows  the  same  idea  as  in  the  isotropic 
analysis,  but  it  uses  Hill's  yield  criterion  and  needs  more  input  as 
far  as  yield  stresses  are  concerned.'  These  yield  stresses  are  in  the 
six  orthotropic  directions  and  are  referred  to  as  Y..  (i  = 1,  3 and  j = 
1,  3).  Before  seeing  the  yield  criterion,  the  following  constants  are 
defined  as: 


1_ 

1 

Y11 

Y2 

Y2 

N 1 
, > 

l 

11 

22 

33 

1 

1_ 

Y22 

Y2 

" Y?, 

" Y2 

22 

11 

33 

1_ 

1_ 

Y33 

= Y2 

" Y2 

Y2 

33 

11 

22 

Then  the  yield  criterion  can  be  written  as: 
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(31) 


f (cr.  .) 
. iJ 


°11  , °22  °33  , °n  ah  °*23 

y2  y2  yZ  y2  y2 

11  22  33  12  13  23 


+ Ylia22a33  + Y22ail°33  + Y33ail022  = 1 


The  total  stress  at  t ..  is: 

n+l 


T T 

a. . = a. . + da. . 

1‘*n+l  1Jn 


n+l 


(32) 


G p 

The  strain  increment  is  divided  into  elastic  (e  ) and  plastic  (e  ) 
strain. 


de.  . = def . + de?. 
ij  13  13 


(33) 


The  flow  rule  is  written  in  the  following  manner: 


den  ’ dX 


11 


V 'll 


de22  = dX 


22 


22 


Y22°33  + Y33°22 


= dXT 


11 


+ Yllg33*Y33°ll,=  dXT 


22 


(34) 


~~33 


= rU 


jl 
V 


♦ Ill?22  * Y22°ll  , . H>T 

2 / " 33 


14 


d£12  = dX 


del3  ’ dX 


de23  = dX 


'12 


12 


13 


13 


23 


13 


= dXT 


12 


= dXT 


13 


dXT 


23 


where  T.^  is  defined  by  Equation  (34).  This  could  be  written  as: 


de?.  = dX  T. . 
ij  13 


(35) 


The  orthotropic  elastic  relation  used  to  evaluate  the  trail  stress  is: 


da.  . = C.  ..  ..  de,  . 
13  ljkl  kl 


(36) 


This  is  substituted  into  Equation  (32)  and  that  result  is  substituted 

T 

into  Equation  (31).  Similar  to  the  previous  analysis  if  f(a. . ) < 1, 
'T  *p  1J  ' ■ 

then  a. . . = a. . , but  if  f(a. . ) > 1 plastic  flow  has  occurred. 

ljn+1  1Jn+l  13 

Inserting  Equation  (33)  into  (36)  and  that  result  into  (32)  gives: 


a.  . = a.  .*  - C.  ...  T,  . dX 
13  i]  i3kl  kl 

Letting  T^*  Equation  (37)  becomes 


(37) 


a. . = a. .T  - T. . dX 
13  iJ  ij 


(38) 


Substituting  Equation  (38)  into  (31)  gives: 

C 


dX 


/ 


(39) 


B + ,/  B - AC 
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where : 


T? 


T2  t"2  'i’2  7f?2  m2 

Y2  Y2  Y2  Y2  v2  v2  11  22  33 

11  *22  *33  *12  Y13  Y23 


(40) 


+ YTT  +YTT 
22  11 1 33  *33  11  22 


B = 


if 11  + T22g22  + T33g33  + flflS + flfll  + T23CT23' 


11 


22 


33 


13 


12 


23 


+ Y 


11 


— T — T 

f T22°33  + T33°22 

Y I 

— T — T 

f Tlia33  + T33ail  1 

1 2 1 

Y22  1 

1 2 ' 

+ Y 


33 


Tlia22T+  T22gllT 


C = f (at j)  - l 


This  dX  is  substituted  into  Equation  (38)  to  give  the  total  stress. 

In  all  cases  this  breaks  down  to  the  isotropic  case  when  using  appro- 
priate yield  stress.  After  programming,  the  same  results  were  produced 
for  an  isotropic  example  as  were  produced  in  the  isotropic  elastic- 
plastic  analysis. 


IV.  VISCO-PLASTIC  MODEL 


Analysis 

In  this  analysis  it  is  assumed  that  the  strain  is  divided  into 
elastic  (dce)  and  visco-plastic  (devP)  strain. 

de  = dee  + de^  (41) 
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Using  Hill's  Flow  Rules,  Equation  (35)  is  modified  to  be: 

de^  = dA  T..  (42) 

It  should  be  noted  that  the  visco-plastic  strain  changes  satisfy  incom- 
pressibility condition 

3 

l de  ? = 0 (43) 

i-1  11 

The  quantities  T-jj  represent  six  independent  quantities  and  they  can  be 
arranged  in  a matrix  form  and  then  can  be  related  to  a stress  matrix  as 
follows: 

Ifl  = TRl  Irrl  (AA^\ 

^ * j j ^ ~r 

where: 
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By  comparing  the  flow  rule  in  equation  (42)  to  an  isotropic  case  it  can 
be  noted  that  the  quantities  T^j  have  the  same  role  as  the  deviatoric 
stresses  and  e^?  strains  as  the  deviatoric  strains.  In  fact,  it  may  be 
noted  that  equation  (42)  reduces  in  the  limit  to  the  isotropic  case. 
Consequently,  following  the  procedure  developed  for  isotropic  Bingham 
material,  the  strain  rate  dependence  is  introduced  by  defining  {Tf} 

{tf}  = {T}  + n (46) 

where  ri  represents  viscous  coefficient,  and  (T)  is  the  quantity  which 
satisfies  the  yield  criterion. 

By  using  Equation  (44)  we  define: 

(aF)  = [B]'1  {TF}  (47) 

If  no  plastic  yield  has  occurred  the  viscoplastic  strain  increment  is 
zero.  So  we  begin  this  analysis  with  a trial  incremental  stress  where: 

{daT}  = [C]  {de}  (48) 


and  [C]  is  the  orthotropic  relationship  between  stress  and  strain. 
Equation  (48)  is  inserted  into  Equation  (32)  and  this  is  inserted  into 
the  yield  criterion  in  Equation  (31).  If  yield  does  not  occur,  the 
trial  stress  is  equal  to  ap.  However,  if  yield  does  occur  then  Equa- 
tions (46)  and  (47)  must  be  used  to  calculate  ap 
Equation  (41)  and  (48)  : 


as  follows:  From 


{doF>  = [C]  {de}  - [C]  {de^} 


(49) 


Then  as  in  Equation  (32): 

{oF}  ={aT}  - [C]  {de^}  (50) 

Multiplying  Equation  (46)  by  [B]  1 and  using  Equation  (47)  in  (46)  one 
gets : 

{aF}  = {a}  + n [B] -1  {evp}  (51) 

From  Equation  (50)  and  (51),  when  solving  for  {a}  one  finds: 

{a}  = {aT}  - [C]  {devp}  - n [B]"1  {e^}  (52) 

R \r  ncinn  Unnof  i nn 

LSJ  UJX11&  l.^uuuj.oh 

{evp}  = 


t +■ 

= 

1 dt  1 


i c?  o o c n "I  \ r e Vi/v.m  f hof  ‘ 
jlu  vciux  j.)  ojiumi  uiac  , 

it {deVP}  - 1 W 


(53) 
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Substituting  (42)  and  (53)  into  (52)  produces: 


{o}  = {oT}  - dX[C]{T}  - [B]_1  {t}  (54) 


Defining  another  variable  (T): 

{T}  = [C]  (T)  (55) 

We  use  the  inverse  of  Equation  (44)  and  Equation  (55)  in  (54)  which 
produces: 


iaj  = 


r T, 

la  i - 


dA 


(r~ 


T •>  \ 


/ r /■  N 


(.OOj 


The  question  in  a dynamic  problem  always  arises  as  to  what  value  of 
stress  is  used  for  the  flow  rule.  In  this  formulation  the  stress  used 
in  the  flow  rule  is  approximated  by  the  trail  stress.  A closer  appro- 
ximation can  be  formed  by  doing  an  iterative  loop  on  this  equation, 
but  little  difference  is  found  in  the  solution  when  this  is  done. 


The  (a)  stress  formed  here  in  Equation  (56)  is  substituted  into 
the  yield  condition  and  dX  is  solved  for. 


This  gives  the  relation: 

AdA2  - 2BdA  + C = 0 


where  letting: 


T.* 

1J 

— T 

= t:.  + 

naT 

At 

A = 

T*2 

-11  ♦ 

T*2 

-22.  + 

T1  *2 

-22.  + 

T*2 

12 

T*2 

+ 13 

T*2 

+ 

Y2 

11 

Y2 

22 

Y2 

33 

Y2 

12 

Y2 

13 

Y2 

23 

+ 

Y T* 
22  11 

T*  + 
33 

Y x* 

11  22 

T* 

33 

+ Y33 

T*  T* 
11  22 

(57) 


(58) 

(59) 
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T*  a T T*  O T T*  a T T*  O T T*  a T T*  a T 
„ _ 11  ll‘  22  22  . 33  331  12  12  . 13  13  . 23  23 

D — + + 11  - + ‘ ""  • + J--  


11 


22 


33 


12 


13 


23 


+ Y 


T T 

T*  a + T*  a 
22  33  33  22 


11 


+ Y 


22 


T T 

T*  a + T*  a 
11  33  33  11 


+ Y 


33 


T T 

T*  a + t*  a 
22  11  11  22 


(60) 


C = f(aj.)  - 1 


(61) 


where  f = 1 is  the  yield  function. 
Therefore,  from  Equation  (57) 


dX 


2B  i /(2B)2  - 4AC 
2A 


(62) 


Since  dX  -*■  0 as  C -►  0 the  minus  sign  must  be  used.  Multiplying  top 
and  bottom  by  B + /b2  - AC  produces: 


dX 


(63) 


This  value  of  dX  is  then  substituted  into  Equation  (56) 
Equations  (51),  (42),  and  (44) 


<aFJ 


1 + 


ndxx 


{o} 


By  using 


(64) 


Equation  (64)  is  used  in  the  modified  computer  program  to  calculate 
the  actual  stress  state  in  finite-elements. 
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Numerical  Example 


The  example  used  to  check  out  the  changes  in  the  program  was  a 
three  layer  laminated  plate4.  The  top  and  bottom  layer  was  1020  steel, 
1,27  cm  and  1.27  cm  thick,  respectively.  The  middle  layer  was  2040 
aluminum  .635  cm  thick.  Values  for  the  viscosity  for  each  material 
were  estimated  to  be: 


298.86  pascal-sec 
18.73  pascal-sec 

The  results  for  this  example  are  given  in  Table  1.  Tabie  1 com- 
pares with  the  vertical  displacement  at  the  center  of  the  plate  of 
both  for  the  original  elastic-plastic  and  the  new  elastic-visco- 
plastic models.  It  can  be  seen  that  for  short  periods  of  time  the 
viscous  effects  are  small,  but  for  later  times  the  effect  is  more  pro- 
nounced. As  expected,  the  effect  of  viscosity  is  to  stiffen  the 
plate. 

The  effect  of  the  viscous  material  properties  do  not  appreciably 
alter  the  dynamic  response  of  steel  and  aluminum  materials.  However, 
this  conclusion  may  not  be  true  for  other  materials,  such  as  for  ex- 
ample, composites,  which  may  exhibit  larger  viscosity. 


V.  CONCLUSIONS 

By  certain  modifications  to  the  TIP  computer  program  it  has  been 
possible  to  increase  the  incremental  time  step  used  in  the  finite 
difference  time  integration.  However,  even  with  these  modifications, 
the  stresses  still  oscillate  in  the  thickness  direction.  It  is  quite 
feasible  that  if  the  oscillations  of  the  stress  can  be  reduced,  fur- 
ther increase  in  time  interval  will  be  possible  with  an  appropriate 
increase  in  numerical  stability. 

The  inclusion  in  the  analysis  and  the  computer  program  of  the 
orthotropic  yield  criterion  and  the  visco-plastic  material  response 
has  been  successfully  accomplished.  The  viscous  effects,  however, 
have  been  found  to  be  small  in  the  case  of  one  numerical  example 
which  uses  aluminum  and  steel  materials. 


4 

Majerus,  J.N.,  and  Knapp,  R.R.,  "Dynamic  Behavior  of  Multi- 
Layered  Plate  Due  to  an  Intense  Impulsive  Load,"  Proceedings  of  the 
Second  International  Conference  on  Mechanical  Behavior  of  Materials, 
Boston,  Mass.,  August,  1976. 
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TOP 


BOTTOM 


Figure  1.  Model  representing  the  predicted  stiffness  in 
the  thickness  direction 
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Figure  2 . Numerical  results  for  time  increment  of 
2.5  us 
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DEFLECTIDN  IN  CENTIMETERS 


Figure  4.  Numerical  results  for  variable  time  increments 
of  1,0,  5.0  and  10.0  us 
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TABLE  1 


Comparison  of  Results  from  the  Elastic-Plastic 
and  Elastic-Visco-Plastic  Models  Using  Variable 
Time  Steps  (see  Figure  4) 


Time 

Vertical  Plate  Displacement 
(cm)  for  Elastic-Plastic 
Model 

Vertical  Plate  Displacement 
(cm)  for  Elastic-Visco- 
Plastic  Model 

10 

-.285747 

-.3158007 

20 

-.666216 

-.6348882 

30 

-.9870059 

-.898324 

40 

-1.277305 

-1.1402441 

50 

-1.5436494 

-1.3339216 

65 

-1.9010985 

-1.5438323 

75 

-2.1077580 

-1.6701541 

85 

-2.2934676 

-1.8299201 

95 

-2,4730964 

-1.9608596 

105 

-2,6440384 

-1.9884898 
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